module operators_mod

  use const_mod
  use perf_mod
  use vert_coord_mod
  use block_mod
  use latlon_parallel_mod
  use latlon_field_types_mod
  use latlon_operators_mod
  use process_mod, only: process_stop
  use formula_mod
  use namelist_mod
  use tracer_mod
  use pgf_mod
  use adv_mod
  use interp_mod
  use filter_mod

  implicit none

  private

  public operators_init
  public operators_prepare
  public calc_mg
  public calc_ph
  public calc_omg
  public calc_dmg
  public calc_t
  public calc_rhod
  public calc_gz
  public calc_mfz
  public calc_div
  public calc_vor
  public calc_coriolis
  public calc_grad_ke
  public calc_grad_mf
  public calc_grad_ptf
  public calc_dmgsdt
  public calc_wedudlev_wedvdlev

  interface operators_prepare
    module procedure operators_prepare_1
    module procedure operators_prepare_2
  end interface operators_prepare

  interface
    subroutine interp_pv_interface(block, dstate, dt, substep)
      import block_type, dstate_type, r8
      type(block_type), intent(inout) :: block
      type(dstate_type), intent(inout) :: dstate
      real(r8), intent(in) :: dt
      integer, intent(in) :: substep
    end subroutine interp_pv_interface
  end interface

  procedure(interp_pv_interface), pointer :: interp_pv => null()

contains

  subroutine operators_init()

    select case (pv_adv_scheme)
    case ('midpoint')
      interp_pv => interp_pv_midpoint
    case ('upwind')
      interp_pv => interp_pv_upwind
    case ('weno')
      interp_pv => interp_pv_weno
    case default
      if (proc%is_root()) call log_error('Invalid pv_scheme ' // trim(pv_adv_scheme) // '!')
    end select

  end subroutine operators_init

  ! First time call before main model loop.

  subroutine operators_prepare_1(blocks, itime, dt)

    type(block_type), intent(inout) :: blocks(:)
    integer, intent(in) :: itime
    real(r8), intent(in) :: dt

    integer iblk

    do iblk = 1, size(blocks)
      if (baroclinic    ) call calc_mg    (blocks(iblk), blocks(iblk)%dstate(itime))
      call calc_dmg                       (blocks(iblk), blocks(iblk)%dstate(itime))
      call tracer_calc_qm                 (blocks(iblk))
      if (baroclinic    ) call calc_ph    (blocks(iblk), blocks(iblk)%dstate(itime))
      if (nonhydrostatic .and. .not. restart) then
        if (sum(blocks(iblk)%dstate(itime)%p%d) == 0) then
          ! Set pressure to hydrostatic pressure in the initial condition.
          blocks(iblk)%dstate(itime)%p    %d = blocks(iblk)%dstate(itime)%ph    %d
          blocks(iblk)%dstate(itime)%p_lev%d = blocks(iblk)%dstate(itime)%ph_lev%d
        end if
      end if
      if (baroclinic    ) call calc_t     (blocks(iblk), blocks(iblk)%dstate(itime))
      call calc_mf                        (blocks(iblk), blocks(iblk)%dstate(itime), dt)
      call calc_ke                        (blocks(iblk), blocks(iblk)%dstate(itime),     total_substeps)
      call calc_pv                        (blocks(iblk), blocks(iblk)%dstate(itime),     total_substeps)
      call interp_pv                      (blocks(iblk), blocks(iblk)%dstate(itime), dt, total_substeps)
      if (baroclinic .and. (hydrostatic .or. init_hydrostatic_gz)) then
        call calc_gz                      (blocks(iblk), blocks(iblk)%dstate(itime))
      end if
      if (baroclinic    ) call calc_rhod  (blocks(iblk), blocks(iblk)%dstate(itime))
    end do

  end subroutine operators_prepare_1

  subroutine operators_prepare_2(block, dstate, dt, pass, substep)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    real(r8), intent(in) :: dt
    integer, intent(in) :: pass
    integer, intent(in) :: substep

    select case (pass)
    case (all_pass)
      call calc_mf                        (block, dstate, dt)
      call calc_ke                        (block, dstate,     substep)
      call calc_pv                        (block, dstate,     substep)
      call interp_pv                      (block, dstate, dt, substep)
      if (baroclinic    ) call calc_t     (block, dstate)
      if (hydrostatic   ) call calc_gz    (block, dstate)
      if (baroclinic    ) call calc_rhod  (block, dstate)
    case (forward_pass)
      call calc_mf                        (block, dstate, dt)
      call calc_ke                        (block, dstate,     substep)
      call calc_pv                        (block, dstate,     substep)
      call interp_pv                      (block, dstate, dt, substep)
    case (backward_pass)
      if (baroclinic    ) call calc_t     (block, dstate)
      if (hydrostatic   ) call calc_gz    (block, dstate)
      if (baroclinic    ) call calc_rhod  (block, dstate)
    end select

  end subroutine operators_prepare_2

  subroutine calc_mg(block, dstate)

    type(block_type), intent(in) :: block
    type(dstate_type), intent(inout) :: dstate

    integer i, j, k, is, ie, js, je

    call perf_start('calc_mg')

    associate (mesh    => block%mesh    , &
               mgs     => dstate%mgs    , & ! in
               mg_lev  => dstate%mg_lev , & ! out
               mg      => dstate%mg     )   ! out
    is = mesh%full_ids - 1
    ie = mesh%full_ide + 1
    js = mesh%full_jds - merge(0, 1, mesh%has_south_pole())
    je = mesh%full_jde + merge(0, 1, mesh%has_north_pole())
    do k = mesh%half_kds, mesh%half_kde
      do j = js, je
        do i = is, ie
          mg_lev%d(i,j,k) = vert_coord_calc_mg_lev(k, mgs%d(i,j), block%static%ref_ps_perb%d(i,j))
        end do
      end do
    end do
    do k = mesh%full_kds, mesh%full_kde
      do j = js, je
        do i = is, ie
          mg%d(i,j,k) = 0.5_r8 * (mg_lev%d(i,j,k) + mg_lev%d(i,j,k+1))
        end do
      end do
    end do
    end associate

    call perf_stop('calc_mg')

  end subroutine calc_mg

  subroutine calc_ph(block, dstate)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate

    integer i, j, k, is, ie, js, je

    call perf_start('calc_ph')

    associate (mesh    => block%mesh          , &
               mg      => dstate%mg           , & ! in
               mg_lev  => dstate%mg_lev       , & ! in
               dmg     => dstate%dmg          , & ! in
               qm      => tracers(block%id)%qm, & ! in
               ph_lev  => dstate%ph_lev       , & ! out
               ph      => dstate%ph           )   ! out
    is = mesh%full_ids - 1
    ie = mesh%full_ide + 1
    js = mesh%full_jds - merge(0, 1, mesh%has_south_pole())
    je = mesh%full_jde + merge(0, 1, mesh%has_north_pole())
    k = mesh%half_kds
    ph_lev%d(:,:,k) = mg_lev%d(:,:,k)
    do k = mesh%half_kds + 1, mesh%half_kde
      do j = js, je
        do i = is, ie
          ph_lev%d(i,j,k) = ph_lev%d(i,j,k-1) + dmg%d(i,j,k-1) * (1 + qm%d(i,j,k-1))
        end do
      end do
    end do
    do k = mesh%full_kds, mesh%full_kde
      do j = js, je
        do i = is, ie
          ph%d(i,j,k) = 0.5_r8 * (ph_lev%d(i,j,k) + ph_lev%d(i,j,k+1))
        end do
      end do
    end do
    end associate

    call perf_stop('calc_ph')

  end subroutine calc_ph

  subroutine calc_omg(block, dstate)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate

    integer i, j, k
    real(r8) sum_dmf(block%mesh%full_ids:block%mesh%full_ide, &
                     block%mesh%full_jds:block%mesh%full_jde)
    real(r8) work   (block%mesh%full_ids:block%mesh%full_ide, &
                     block%mesh%full_kds:block%mesh%full_kde)
    real(r8) pole   (block%mesh%full_kds:block%mesh%full_kde)

    call perf_start('calc_omg')

    associate (mesh  => block%mesh   , &
               ph    => dstate%ph    , & ! in
               u_lon => dstate%u_lon , & ! in
               v_lat => dstate%v_lat , & ! in
               dmf   => block%aux%dmf, & ! in
               div   => block%aux%div, & ! in
               omg   => block%aux%omg)   ! out
    sum_dmf = 0
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
        do i = mesh%full_ids, mesh%full_ide
          sum_dmf(i,j) = sum_dmf(i,j) + dmf%d(i,j,k)
          omg%d(i,j,k) = 0.5_r8 * ((                                              &
            u_lon%d(i  ,j,k) * (ph%d(i,j,k) + ph%d(i+1,j,k)) -                    &
            u_lon%d(i-1,j,k) * (ph%d(i,j,k) + ph%d(i-1,j,k))                      &
          ) * mesh%le_lon(j) + (                                                  &
            v_lat%d(i,j  ,k) * (ph%d(i,j,k) + ph%d(i,j+1,k)) * mesh%le_lat(j  ) - &
            v_lat%d(i,j-1,k) * (ph%d(i,j,k) + ph%d(i,j-1,k)) * mesh%le_lat(j-1)   &
          )) / mesh%area_cell(j) - ph%d(i,j,k) * div%d(i,j,k) - sum_dmf(i,j)
        end do
      end do
    end do
    if (mesh%has_south_pole()) then
      j = mesh%full_jds
      do k = mesh%full_kds, mesh%full_kde
        do i = mesh%full_ids, mesh%full_ide
          work(i,k) = 0.5_r8 * v_lat%d(i,j,k) * (ph%d(i,j,k) + ph%d(i,j+1,k)) * mesh%le_lat(j)
        end do
      end do
      call zonal_sum(proc%zonal_circle, work, pole)
      pole = pole / mesh%area_pole_cap
      do k = mesh%full_kds, mesh%full_kde
        i = mesh%full_ids
        sum_dmf(i,j) = sum_dmf(i,j) + dmf%d(i,j,k)
        omg%d(:,j,k) = pole(k) - ph%d(i,j,k) * div%d(i,j,k) - sum_dmf(i,j)
      end do
    end if
    if (mesh%has_north_pole()) then
      j = mesh%full_jde
      do k = mesh%full_kds, mesh%full_kde
        do i = mesh%full_ids, mesh%full_ide
          work(i,k) = 0.5_r8 * v_lat%d(i,j-1,k) * (ph%d(i,j,k) + ph%d(i,j-1,k)) * mesh%le_lat(j-1)
        end do
      end do
      call zonal_sum(proc%zonal_circle, work, pole)
      pole = -pole / mesh%area_pole_cap
      do k = mesh%full_kds, mesh%full_kde
        i = mesh%full_ids
        sum_dmf(i,j) = sum_dmf(i,j) + dmf%d(i,j,k)
        omg%d(:,j,k) = pole(k) - ph%d(i,j,k) * div%d(i,j,k) - sum_dmf(i,j)
      end do
    end if
    end associate

    call perf_stop('calc_omg')

  end subroutine calc_omg

  subroutine calc_t(block, dstate)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate

    integer i, j, k

    call wait_halo(dstate%pt)

    call perf_start('calc_t')

    associate (mesh => block%mesh         , &
               pt   => dstate%pt          , & ! in
               p    => dstate%p           , & ! in
               q    => tracers(block%id)%q, & ! in
               t    => block%aux%t        , & ! out
               tv   => block%aux%tv       )   ! out
    if (idx_qv > 0) then
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds, mesh%full_jde + merge(0, 1, mesh%has_north_pole())
          do i = mesh%full_ids, mesh%full_ide + 1
            t%d(i,j,k) = temperature(pt%d(i,j,k), p%d(i,j,k), q%d(i,j,k,idx_qv))
            tv%d(i,j,k) = virtual_temperature_from_modified_potential_temperature(pt%d(i,j,k), p%d(i,j,k)**rd_o_cpd, q%d(i,j,k,idx_qv))
          end do
        end do
      end do
    else
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds, mesh%full_jde + merge(0, 1, mesh%has_north_pole())
          do i = mesh%full_ids, mesh%full_ide + 1
            t%d(i,j,k) = temperature(pt%d(i,j,k), p%d(i,j,k), 0.0_r8)
            tv%d(i,j,k) = t%d(i,j,k)
          end do
        end do
      end do
    end if
    end associate

    call perf_stop('calc_t')

  end subroutine calc_t

  subroutine calc_rhod(block, dstate)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate

    integer i, j, k

    call wait_halo(dstate%gz_lev)

    call perf_start('calc_rhod')

    associate (mesh   => block%mesh    , &
               gz_lev => dstate%gz_lev , & ! in
               dmg    => dstate%dmg    , & ! in
               rhod   => block%aux%rhod)   ! out
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds, mesh%full_jde + merge(0, 1, mesh%has_north_pole())
        do i = mesh%full_ids, mesh%full_ide + 1
          rhod%d(i,j,k) = dmg%d(i,j,k) / (gz_lev%d(i,j,k) - gz_lev%d(i,j,k+1))
        end do
      end do
    end do
    end associate

    call perf_stop('calc_rhod')

  end subroutine calc_rhod

  subroutine calc_mfz(block, dstate, dtend)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(in) :: dstate
    type(dtend_type), intent(in) :: dtend

    integer i, j, k

    call perf_start('calc_mfz')

    associate (mesh        => block%mesh       , &
               dmf         => block%aux%dmf    , & ! in
               dmgsdt      => dtend%dmgsdt     , & ! in
               mfz_lev     => block%aux%mfz_lev)   ! out
    do k = mesh%half_kds + 1, mesh%half_kde - 1
      do j = mesh%full_jds, mesh%full_jde
        do i = mesh%full_ids, mesh%full_ide
          mfz_lev%d(i,j,k) = -vert_coord_calc_dmgdt_lev(k, dmgsdt%d(i,j)) - sum(dmf%d(i,j,1:k-1))
        end do
      end do
    end do
    call fill_halo(mfz_lev, west_halo=.false., south_halo=.false., async=.true.)
    end associate

    call perf_stop('calc_mfz')

  end subroutine calc_mfz

  subroutine calc_ke(block, dstate, substep)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    integer, intent(in) :: substep

    integer i, j, k
    real(r8) ke_vtx(4)
    real(r8) work(block%mesh%full_ids:block%mesh%full_ide,block%mesh%full_nlev)
    real(r8) pole(block%mesh%full_nlev)

    call perf_start('calc_ke')

    associate (mesh => block%mesh  , &
               u    => dstate%u_lon, & ! in
               v    => dstate%v_lat, & ! in
               ke   => block%aux%ke)   ! out
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole + merge(0, 1, mesh%has_north_pole())
        do i = mesh%full_ids, mesh%full_ide + 1
          ke%d(i,j,k) = 0.25_r8 * (             &
            u%d(i-1,j  ,k)**2 + u%d(i,j,k)**2 + &
            v%d(i  ,j-1,k)**2 + v%d(i,j,k)**2   &
          )
        end do
      end do
    end do
    if (mesh%has_south_pole()) then
      j = mesh%full_jds
      do k = mesh%full_kds, mesh%full_kde
        do i = mesh%full_ids, mesh%full_ide
          work(i,k) = v%d(i,j,k)**2
        end do
      end do
      call zonal_sum(proc%zonal_circle, work, pole)
      pole = pole / global_mesh%full_nlon
      do k = mesh%full_kds, mesh%full_kde
        do i = mesh%full_ids, mesh%full_ide
          ke%d(i,j,k) = pole(k)
        end do
      end do
    end if
    if (mesh%has_north_pole()) then
      j = mesh%full_jde
      do k = mesh%full_kds, mesh%full_kde
        do i = mesh%full_ids, mesh%full_ide
          work(i,k) = v%d(i,j-1,k)**2
        end do
      end do
      call zonal_sum(proc%zonal_circle, work, pole)
      pole = pole / global_mesh%full_nlon
      do k = mesh%full_kds, mesh%full_kde
        do i = mesh%full_ids, mesh%full_ide
          ke%d(i,j,k) = pole(k)
        end do
      end do
    end if

    if (save_dyn_calc .and. substep < total_substeps) then
      call perf_stop('calc_ke')
      return
    end if

    if (ke_scheme == 2) then
      !
      !      ________u_________________u________
      !     |     i-1,j+1     |       i,j+1     |
      !     |                 |                 |
      !     |                 |                 |
      !     |        1        |        4        |
      !     v        o--------v--------o        v
      !  i-1,j    i-1,j      i,j      i,j    i+1,j
      !     |        |        |        |        |
      !     |        |        |        |        |
      !     |________u________|________u________|
      !     |     i-1,j      i,j      i,j       |
      !     |        |        |        |        |
      !     |        |        |        |        |
      !     |        |        |        |        |
      !     v        o--------v--------o        v
      !  i-1,j-1  i-1,j-1    i,j-1    i,j-1  i+1,j-1
      !     |        2        |        3        |
      !     |                 |                 |
      !     |________u________|________u________|
      !           i-1,j-1             i,j-1
      !
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole + merge(0, 1, mesh%has_north_pole())
          do i = mesh%full_ids, mesh%full_ide + 1
            ke_vtx(1) = 0.25_r8 * (                   &
              v%d(i-1,j  ,k)**2 + v%d(i  ,j  ,k)**2 + &
              u%d(i-1,j  ,k)**2 + u%d(i-1,j+1,k)**2   &
            )
            ke_vtx(2) = 0.25_r8 * (                   &
              v%d(i-1,j-1,k)**2 + v%d(i  ,j-1,k)**2 + &
              u%d(i-1,j-1,k)**2 + u%d(i-1,j  ,k)**2   &
            )
            ke_vtx(3) = 0.25_r8 * (                   &
              v%d(i  ,j-1,k)**2 + v%d(i+1,j-1,k)**2 + &
              u%d(i  ,j-1,k)**2 + u%d(i  ,j  ,k)**2   &
            )
            ke_vtx(4) = 0.25_r8 * (                   &
              v%d(i  ,j  ,k)**2 + v%d(i+1,j  ,k)**2 + &
              u%d(i  ,j  ,k)**2 + u%d(i  ,j+1,k)**2   &
            )
            ke%d(i,j,k) = (1.0_r8 - ke_cell_wgt) * 0.25_r8 * ( &
              ke_vtx(1) + ke_vtx(4) + ke_vtx(2) + ke_vtx(3)    &
            ) + ke_cell_wgt * ke%d(i,j,k)
          end do
        end do
      end do
    end if
    end associate

    call perf_stop('calc_ke')

  end subroutine calc_ke

  subroutine calc_div(block, dstate)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate

    integer i, j, k

    call wait_halo(dstate%u_lon)
    call wait_halo(dstate%v_lat)

    call perf_start('calc_div')

    associate (mesh  => block%mesh        , &
               u_lon => dstate%u_lon      , & ! in
               v_lat => dstate%v_lat      , & ! in
               div   => block%aux%div     , & ! out
               divx  => block%aux%g_3d_lon, & ! working array
               divy  => block%aux%g_3d_lat, & ! working array
               div2  => block%aux%div2    )   ! out
    call div_operator(u_lon, v_lat, div, with_halo=.true.)
    if (div_damp_order == 4) then
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
          do i = mesh%half_ids, mesh%half_ide + 1
            divx%d(i,j,k) = (div%d(i+1,j,k) - div%d(i,j,k)) / mesh%de_lon(j)
          end do
        end do
      end do
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%half_jds, mesh%half_jde + merge(0, 1, mesh%has_north_pole())
          do i = mesh%full_ids, mesh%full_ide
            divy%d(i,j,k) = (div%d(i,j+1,k) - div%d(i,j,k)) / mesh%de_lat(j)
          end do
        end do
      end do
      call div_operator(divx, divy, div2)
      call filter_run(block%big_filter, div2)
      call fill_halo(div2, west_halo=.false., south_halo=.false.)
    else
      call filter_run(block%small_filter, div)
      call fill_halo(div, west_halo=.false., south_halo=.false.)
    end if
    end associate

    call perf_stop('calc_div')

  end subroutine calc_div

  subroutine calc_gz(block, dstate)

    type(block_type), intent(in) :: block
    type(dstate_type), intent(inout) :: dstate

    integer i, j, k

    call perf_start('calc_gz')

    associate (mesh   => block%mesh      , &
               tv     => block%aux%tv    , & ! in
               ph_lev => dstate%ph_lev   , & ! in
               gz_lev => dstate%gz_lev   , & ! out
               gz     => dstate%gz       )   ! out
    do k = mesh%half_kde - 1, mesh%half_kds, -1
      do j = mesh%full_jds, mesh%full_jde + merge(0, 1, mesh%has_north_pole())
        do i = mesh%full_ids, mesh%full_ide + 1
          gz_lev%d(i,j,k) = gz_lev%d(i,j,k+1) + rd * tv%d(i,j,k) * log(ph_lev%d(i,j,k+1) / ph_lev%d(i,j,k))
        end do
      end do
    end do
    ! For output
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds, mesh%full_jde + merge(0, 1, mesh%has_north_pole())
        do i = mesh%full_ids, mesh%full_ide + 1
          gz%d(i,j,k) = 0.5_r8 * (gz_lev%d(i,j,k) + gz_lev%d(i,j,k+1))
        end do
      end do
    end do
    end associate

    call perf_stop('calc_gz')

  end subroutine calc_gz

  subroutine calc_dmg(block, dstate)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate

    integer i, j, k, l, is, ie, js, je

    call perf_start('calc_dmg')

    associate (mesh    => block%mesh       , &
               mg      => dstate%mg        , & ! in
               mg_lev  => dstate%mg_lev    , & ! in
               gz      => dstate%gz        , & ! in
               gzs     => block%static%gzs , & ! in
               dmg     => dstate%dmg       , & ! out
               dmg_lon => block%aux%dmg_lon, & ! out
               dmg_lat => block%aux%dmg_lat, & ! out
               dmg_lev => dstate%dmg_lev   , & ! out
               dmg_vtx => block%aux%dmg_vtx)   ! out
    is = mesh%full_ids - 1
    ie = mesh%full_ide + 1
    js = mesh%full_jds - merge(0, 1, mesh%has_south_pole())
    je = mesh%full_jde + merge(0, 1, mesh%has_north_pole())
    if (baroclinic .or. advection) then
      do k = mesh%full_kds, mesh%full_kde
        do j = js, je
          do i = is, ie
            dmg%d(i,j,k) = mg_lev%d(i,j,k+1) - mg_lev%d(i,j,k)
            if (dmg%d(i,j,k) <= 0) then
              do l = mesh%half_kds, mesh%half_kde
                print *, l, mg_lev%d(i,j,l)
              end do
              print *, 'mgs(i,j) =', dstate%mgs%d(i,j)
              print *, mesh%full_lon_deg(i), '(', to_str(i), ')', mesh%full_lat_deg(j), '(', to_str(j), ')', k
              call log_warning('The dry-air weight levels are not monotonic!', __FILE__, __LINE__)
              call process_stop(1)
            end if
          end do
        end do
      end do
      do k = mesh%half_kds + 1, mesh%half_kde - 1
        do j = js, je
          do i = is, ie
            dmg_lev%d(i,j,k) = mg%d(i,j,k) - mg%d(i,j,k-1)
          end do
        end do
      end do
      ! Top boundary
      k = mesh%half_kds
      do j = js, je
        do i = is, ie
          dmg_lev%d(i,j,k) = mg%d(i,j,k) - mg_lev%d(i,j,k)
        end do
      end do
      ! Bottom boundary
      k = mesh%half_kde
      do j = js, je
        do i = is, ie
          dmg_lev%d(i,j,k) = mg_lev%d(i,j,k) - mg%d(i,j,k-1)
        end do
      end do
      call fill_halo(dmg_lev)
    else
      do j = js, je
        do i = is, ie
          dmg%d(i,j,1) = (gz%d(i,j,1) - gzs%d(i,j)) / g
        end do
      end do
    end if

    call average_run(dmg, dmg_lon); call fill_halo(dmg_lon, async=.true.)
    call average_run(dmg, dmg_lat); call fill_halo(dmg_lat, async=.true.)
    call interp_run (dmg, dmg_vtx)
    end associate

    call perf_stop('calc_dmg')

  end subroutine calc_dmg

  subroutine calc_mf(block, dstate, dt)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    real(r8), intent(in) :: dt

    integer i, j, k

    call wait_halo(dstate%u_lon)
    call wait_halo(dstate%v_lat)
    call wait_halo(block%aux%dmg_lon)
    call wait_halo(block%aux%dmg_lat)

    call perf_start('calc_mf')

    associate (mesh    => block%mesh       , &
               dmg     => dstate%dmg       , & ! in
               dmg_lon => block%aux%dmg_lon, & ! in
               dmg_lat => block%aux%dmg_lat, & ! in
               u_lon   => dstate%u_lon     , & ! in
               v_lat   => dstate%v_lat     , & ! in
               u_lat   => block%aux%u_lat  , & ! out
               v_lon   => block%aux%v_lon  , & ! out
               mfx_lon => block%aux%mfx_lon, & ! out
               mfy_lat => block%aux%mfy_lat, & ! out
               mfy_lon => block%aux%mfy_lon, & ! out
               mfx_lat => block%aux%mfx_lat)   ! out
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole + merge(0, 1, mesh%has_north_pole())
        do i = mesh%half_ids - 1, mesh%half_ide
          mfx_lon%d(i,j,k) = dmg_lon%d(i,j,k) * u_lon%d(i,j,k)
        end do
      end do
    end do
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%half_jds - merge(0, 1, mesh%has_south_pole()), mesh%half_jde
        do i = mesh%full_ids, mesh%full_ide + 1
          mfy_lat%d(i,j,k) = dmg_lat%d(i,j,k) * v_lat%d(i,j,k)
        end do
      end do
    end do

    call interp_run(mfx_lon, mfx_lat)
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%half_jds, mesh%half_jde
        do i = mesh%full_ids, mesh%full_ide
          u_lat%d(i,j,k) = mfx_lat%d(i,j,k) / dmg_lat%d(i,j,k)
        end do
      end do
    end do

    call interp_run(mfy_lat, mfy_lon)
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
        do i = mesh%half_ids, mesh%half_ide
          v_lon%d(i,j,k) = mfy_lon%d(i,j,k) / dmg_lon%d(i,j,k)
        end do
      end do
    end do
    end associate

    call perf_stop('calc_mf')

  end subroutine calc_mf

  subroutine calc_vor(block, dstate, with_halo)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    logical, intent(in), optional :: with_halo

    integer i, j, k
    real(r8) work(block%mesh%half_ids:block%mesh%half_ide,block%mesh%full_nlev)
    real(r8) pole(block%mesh%full_nlev)

    call wait_halo(dstate%u_lon)
    call wait_halo(dstate%v_lat)

    call perf_start('calc_vor')

    associate (mesh  => block%mesh     , &
               u_lon => dstate%u_lon   , & ! in
               v_lat => dstate%v_lat   , & ! in
               u_lat => block%aux%u_lat, & ! in
               vor   => block%aux%vor  )   ! out
    call curl_operator(u_lon, v_lat, vor, with_halo)
    if (pv_pole_stokes) then
      ! Special treatment of vorticity around Poles
      if (mesh%has_south_pole()) then
        j = mesh%half_jds
        do k = mesh%full_kds, mesh%full_kde
          do i = mesh%half_ids, mesh%half_ide
            work(i,k) = u_lat%d(i,j,k)
          end do
        end do
        call zonal_sum(proc%zonal_circle, work, pole)
        pole = -pole * mesh%le_lat(j) / global_mesh%area_pole_cap
        do k = mesh%full_kds, mesh%full_kde
          do i = mesh%half_ids - 1, mesh%half_ide
            vor%d(i,j,k) = pole(k)
          end do
        end do
      end if
      if (mesh%has_north_pole()) then
        j = mesh%half_jde
        do k = mesh%full_kds, mesh%full_kde
          do i = mesh%half_ids, mesh%half_ide
            work(i,k) = u_lat%d(i,j,k)
          end do
        end do
        call zonal_sum(proc%zonal_circle, work, pole)
        pole = pole * mesh%le_lat(j) / global_mesh%area_pole_cap
        do k = mesh%full_kds, mesh%full_kde
          do i = mesh%half_ids - 1, mesh%half_ide
            vor%d(i,j,k) = pole(k)
          end do
        end do
      end if
    end if
    end associate

    call perf_stop('calc_vor')

  end subroutine calc_vor

  subroutine calc_pv(block, dstate, substep)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    integer, intent(in) :: substep

    integer i, j, k

    call perf_start('calc_pv')

    associate (mesh    => block%mesh       , &
               dmg_vtx => block%aux%dmg_vtx, & ! in
               vor     => block%aux%vor    , & ! in
               pv      => block%aux%pv     )   ! out
    call calc_vor(block, dstate)
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%half_jds, mesh%half_jde
        do i = mesh%half_ids, mesh%half_ide
          pv%d(i,j,k) = (vor%d(i,j,k) + mesh%f_lat(j)) / dmg_vtx%d(i,j,k)
        end do
      end do
    end do
    call fill_halo(pv)
    end associate

    call perf_stop('calc_pv')

  end subroutine calc_pv

  subroutine interp_pv_midpoint(block, dstate, dt, substep)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    real(r8), intent(in) :: dt
    integer, intent(in) :: substep

    integer i, j, k

    call perf_start('interp_pv_midpoint')

    associate (mesh   => block%mesh      , &
               pv     => block%aux%pv    , & ! in
               pv_lon => block%aux%pv_lon, & ! out
               pv_lat => block%aux%pv_lat)   ! out
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%half_jds, mesh%half_jde
        do i = mesh%full_ids, mesh%full_ide
          pv_lat%d(i,j,k) = 0.5_r8 * (pv%d(i-1,j,k) + pv%d(i,j,k))
        end do
      end do
    end do
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
        do i = mesh%half_ids, mesh%half_ide
          pv_lon%d(i,j,k) = 0.5_r8 * (pv%d(i,j,k) + pv%d(i,j-1,k))
        end do
      end do
    end do
    call fill_halo(pv_lon, east_halo=.false., south_halo=.false., async=.true.)
    call fill_halo(pv_lat, west_halo=.false., north_halo=.false., async=.true.)
    end associate

    call perf_stop('interp_pv_midpoint')

  end subroutine interp_pv_midpoint

  subroutine interp_pv_upwind(block, dstate, dt, substep)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    real(r8), intent(in) :: dt
    integer, intent(in) :: substep

    real(r8) b
    integer i, j, k

    call perf_start('interp_pv_upwind')

    associate (mesh   => block%mesh      , &
               un     => dstate%u_lon    , & ! in
               vn     => dstate%v_lat    , & ! in
               ut     => block%aux%u_lat , & ! in
               vt     => block%aux%v_lon , & ! in
               pv     => block%aux%pv    , & ! in
               pv_lon => block%aux%pv_lon, & ! out
               pv_lat => block%aux%pv_lat)   ! out
    select case (upwind_order_pv)
    case (1)
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
          do i = mesh%half_ids, mesh%half_ide
            b = abs(vt%d(i,j,k)) / (sqrt(un%d(i,j,k)**2 + vt%d(i,j,k)**2) + eps)
            pv_lon%d(i,j,k) = b * upwind1(sign(1.0_r8, vt%d(i,j,k)), upwind_wgt_pv, pv%d(i,j-1:j,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i,j-1,k) + pv%d(i,j,k))
          end do
        end do
        do j = mesh%half_jds, mesh%half_jde
          do i = mesh%full_ids, mesh%full_ide
            b = abs(ut%d(i,j,k)) / (sqrt(ut%d(i,j,k)**2 + vn%d(i,j,k)**2) + eps)
            pv_lat%d(i,j,k) = b * upwind1(sign(1.0_r8, ut%d(i,j,k)), upwind_wgt_pv, pv%d(i-1:i,j,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i-1,j,k) + pv%d(i,j,k))
          end do
        end do
      end do
    case (3)
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
          do i = mesh%half_ids, mesh%half_ide
            b = abs(vt%d(i,j,k)) / (sqrt(un%d(i,j,k)**2 + vt%d(i,j,k)**2) + eps)
            pv_lon%d(i,j,k) = b * upwind3(sign(1.0_r8, vt%d(i,j,k)), upwind_wgt_pv, pv%d(i,j-2:j+1,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i,j-1,k) + pv%d(i,j,k))
          end do
        end do
        do j = mesh%half_jds, mesh%half_jde
          do i = mesh%full_ids, mesh%full_ide
            b  = abs(ut%d(i,j,k)) / (sqrt(ut%d(i,j,k)**2 + vn%d(i,j,k)**2) + eps)
            pv_lat%d(i,j,k) = b * upwind3(sign(1.0_r8, ut%d(i,j,k)), upwind_wgt_pv, pv%d(i-2:i+1,j,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i-1,j,k) + pv%d(i,j,k))
          end do
        end do
      end do
    case (5)
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
          do i = mesh%half_ids, mesh%half_ide
            b = abs(vt%d(i,j,k)) / (sqrt(un%d(i,j,k)**2 + vt%d(i,j,k)**2) + eps)
            pv_lon%d(i,j,k) = b * upwind5(sign(1.0_r8, vt%d(i,j,k)), upwind_wgt_pv, pv%d(i,j-3:j+2,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i,j-1,k) + pv%d(i,j,k))
          end do
        end do
        do j = mesh%half_jds, mesh%half_jde
          do i = mesh%full_ids, mesh%full_ide
            b = abs(ut%d(i,j,k)) / (sqrt(ut%d(i,j,k)**2 + vn%d(i,j,k)**2) + eps)
            pv_lat%d(i,j,k) = b * upwind5(sign(1.0_r8, ut%d(i,j,k)), upwind_wgt_pv, pv%d(i-3:i+2,j,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i-1,j,k) + pv%d(i,j,k))
          end do
        end do
      end do
    end select
    if (substep == total_substeps .or. .not. save_dyn_calc) then
      call fill_halo(pv_lon, east_halo=.false., south_halo=.false., async=.true.)
      call fill_halo(pv_lat, west_halo=.false., north_halo=.false., async=.true.)
    end if
    end associate

    call perf_stop('interp_pv_upwind')

  end subroutine interp_pv_upwind

  subroutine interp_pv_weno(block, dstate, dt, substep)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    real(r8), intent(in) :: dt
    integer, intent(in) :: substep

    real(r8) b
    integer i, j, k

    call perf_start('interp_pv_weno')

    associate (mesh   => block%mesh      , &
               un     => dstate%u_lon    , & ! in
               vn     => dstate%v_lat    , & ! in
               ut     => block%aux%u_lat , & ! in
               vt     => block%aux%v_lon , & ! in
               pv     => block%aux%pv    , & ! in
               pv_lon => block%aux%pv_lon, & ! out
               pv_lat => block%aux%pv_lat)   ! out
    select case (weno_order_pv)
    case (1)
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
          do i = mesh%half_ids, mesh%half_ide
            b = abs(vt%d(i,j,k)) / (sqrt(un%d(i,j,k)**2 + vt%d(i,j,k)**2) + eps)
            pv_lon%d(i,j,k) = b * upwind1(sign(1.0_r8, vt%d(i,j,k)), upwind_wgt_pv, pv%d(i,j-1:j,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i,j-1,k) + pv%d(i,j,k))
          end do
        end do
        do j = mesh%half_jds, mesh%half_jde
          do i = mesh%full_ids, mesh%full_ide
            b = abs(ut%d(i,j,k)) / (sqrt(ut%d(i,j,k)**2 + vn%d(i,j,k)**2) + eps)
            pv_lat%d(i,j,k) = b * upwind1(sign(1.0_r8, ut%d(i,j,k)), upwind_wgt_pv, pv%d(i-1:i,j,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i-1,j,k) + pv%d(i,j,k))
          end do
        end do
      end do
    case (3)
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
          do i = mesh%half_ids, mesh%half_ide
            b = abs(vt%d(i,j,k)) / (sqrt(un%d(i,j,k)**2 + vt%d(i,j,k)**2) + eps)
            pv_lon%d(i,j,k) = b * weno3(sign(1.0_r8, vt%d(i,j,k)), pv%d(i,j-2:j+1,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i,j-1,k) + pv%d(i,j,k))
          end do
        end do
        do j = mesh%half_jds, mesh%half_jde
          do i = mesh%full_ids, mesh%full_ide
            b = abs(ut%d(i,j,k)) / (sqrt(ut%d(i,j,k)**2 + vn%d(i,j,k)**2) + eps)
            pv_lat%d(i,j,k) = b * weno3(sign(1.0_r8, ut%d(i,j,k)), pv%d(i-2:i+1,j,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i-1,j,k) + pv%d(i,j,k))
          end do
        end do
      end do
    case (5)
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
          do i = mesh%half_ids, mesh%half_ide
            b = abs(vt%d(i,j,k)) / (sqrt(un%d(i,j,k)**2 + vt%d(i,j,k)**2) + eps)
            pv_lon%d(i,j,k) = b * weno5(sign(1.0_r8, vt%d(i,j,k)), pv%d(i,j-3:j+2,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i,j-1,k) + pv%d(i,j,k))
          end do
        end do
        do j = mesh%half_jds, mesh%half_jde
          do i = mesh%full_ids, mesh%full_ide
            b = abs(ut%d(i,j,k)) / (sqrt(ut%d(i,j,k)**2 + vn%d(i,j,k)**2) + eps)
            pv_lat%d(i,j,k) = b * weno5(sign(1.0_r8, ut%d(i,j,k)), pv%d(i-3:i+2,j,k)) + &
                              (1 - b) * 0.5_r8 * (pv%d(i-1,j,k) + pv%d(i,j,k))
          end do
        end do
      end do
    end select
    if (substep == total_substeps .or. .not. save_dyn_calc) then
      call fill_halo(pv_lon, east_halo=.false., south_halo=.false., async=.true.)
      call fill_halo(pv_lat, west_halo=.false., north_halo=.false., async=.true.)
    end if
    end associate

    call perf_stop('interp_pv_weno')

  end subroutine interp_pv_weno

  subroutine calc_coriolis(block, dstate, dtend, dt, substep)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    type(dtend_type), intent(inout) :: dtend
    real(r8), intent(in) :: dt
    integer, intent(in) :: substep

    real(r8) tmp
    integer i, j, k

    call wait_halo(block%aux%pv_lon)
    call wait_halo(block%aux%pv_lat)

    call perf_start('calc_coriolis')

    associate (mesh    => block%mesh       , &
               mfx_lon => block%aux%mfx_lon, & ! in
               mfy_lat => block%aux%mfy_lat, & ! in
               mfy_lon => block%aux%mfy_lon, & ! in
               mfx_lat => block%aux%mfx_lat, & ! in
               pv_lon  => block%aux%pv_lon , & ! in
               pv_lat  => block%aux%pv_lat , & ! in
               dudt    => dtend%dudt       , & ! out
               dvdt    => dtend%dvdt       )   ! out
    select case (coriolis_scheme)
    case (1)
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%half_jds, mesh%half_jde
          do i = mesh%full_ids, mesh%full_ide
            tmp = - (                                                            &
              mesh%tg_wgt_lat(1,j) * (                                           &
                mfx_lon%d(i-1,j  ,k) * (pv_lat%d(i,j,k) + pv_lon%d(i-1,j  ,k)) + &
                mfx_lon%d(i  ,j  ,k) * (pv_lat%d(i,j,k) + pv_lon%d(i  ,j  ,k))   &
              ) +                                                                &
              mesh%tg_wgt_lat(2,j) * (                                           &
                mfx_lon%d(i-1,j+1,k) * (pv_lat%d(i,j,k) + pv_lon%d(i-1,j+1,k)) + &
                mfx_lon%d(i  ,j+1,k) * (pv_lat%d(i,j,k) + pv_lon%d(i  ,j+1,k))   &
              )                                                                  &
            ) * 0.5_r8
            dvdt%d(i,j,k) = dvdt%d(i,j,k) + tmp
#ifdef OUTPUT_H1_DTEND
            dtend%dvdt_coriolis%d(i,j,k) = tmp
#endif
          end do
        end do
      end do
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
          do i = mesh%half_ids, mesh%half_ide
            tmp = (                                                              &
              mesh%tg_wgt_lon(1,j) * (                                           &
                mfy_lat%d(i  ,j-1,k) * (pv_lon%d(i,j,k) + pv_lat%d(i  ,j-1,k)) + &
                mfy_lat%d(i+1,j-1,k) * (pv_lon%d(i,j,k) + pv_lat%d(i+1,j-1,k))   &
              ) +                                                                &
              mesh%tg_wgt_lon(2,j) * (                                           &
                mfy_lat%d(i  ,j  ,k) * (pv_lon%d(i,j,k) + pv_lat%d(i  ,j  ,k)) + &
                mfy_lat%d(i+1,j  ,k) * (pv_lon%d(i,j,k) + pv_lat%d(i+1,j  ,k))   &
              )                                                                  &
            ) * 0.5_r8
            dudt%d(i,j,k) = dudt%d(i,j,k) + tmp
#ifdef OUTPUT_H1_DTEND
            dtend%dudt_coriolis%d(i,j,k) = tmp
#endif
          end do
        end do
      end do
    case (2)
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%half_jds, mesh%half_jde
          do i = mesh%full_ids, mesh%full_ide
            dvdt%d(i,j,k) = dvdt%d(i,j,k) - mfx_lat%d(i,j,k) * pv_lat%d(i,j,k)
          end do
        end do
      end do
      do k = mesh%full_kds, mesh%full_kde
        do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
          do i = mesh%half_ids, mesh%half_ide
            dudt%d(i,j,k) = dudt%d(i,j,k) + mfy_lon%d(i,j,k) * pv_lon%d(i,j,k)
          end do
        end do
      end do
    end select
    end associate

    call perf_stop('calc_coriolis')

  end subroutine calc_coriolis

  subroutine calc_grad_ke(block, dstate, dtend, dt)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    type(dtend_type), intent(inout) :: dtend
    real(r8), intent(in) :: dt

    real(r8) tmp
    integer i, j, k

    call perf_start('calc_grad_ke')

    associate (mesh => block%mesh  , &
               ke   => block%aux%ke, & ! in
               dudt => dtend%dudt  , & ! out
               dvdt => dtend%dvdt  )   ! out
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
        do i = mesh%half_ids, mesh%half_ide
          tmp = -(ke%d(i+1,j,k) - ke%d(i,j,k)) / mesh%de_lon(j)
          dudt%d(i,j,k) = dudt%d(i,j,k) + tmp
#ifdef OUTPUT_H1_DTEND
          dtend%dudt_dkedx%d(i,j,k) = tmp
#endif
        end do
      end do
    end do
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%half_jds, mesh%half_jde
        do i = mesh%full_ids, mesh%full_ide
          tmp = -(ke%d(i,j+1,k) - ke%d(i,j,k)) / mesh%de_lat(j)
          dvdt%d(i,j,k) = dvdt%d(i,j,k) + tmp
#ifdef OUTPUT_H1_DTEND
          dtend%dvdt_dkedy%d(i,j,k) = tmp
#endif
        end do
      end do
    end do
    end associate

    call perf_stop('calc_grad_ke')

  end subroutine calc_grad_ke

  subroutine calc_grad_mf(block, dstate)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(in) :: dstate

    call perf_start('calc_grad_mf')

    associate (mesh    => block%mesh       , &
               mfx_lon => block%aux%mfx_lon, & ! in
               mfy_lat => block%aux%mfy_lat, & ! in
               dmf     => block%aux%dmf    )   ! out
    call div_operator(mfx_lon, mfy_lat, dmf)
    end associate

    call perf_stop('calc_grad_mf')

  end subroutine calc_grad_mf

  subroutine calc_grad_ptf(block, dstate, dtend, dt)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(inout) :: dstate
    type(dtend_type ), intent(inout) :: dtend
    real(r8), intent(in) :: dt

    integer i, j, k

    call wait_halo(dstate%pt)

    call perf_start('calc_grad_ptf')

    associate (mesh    => block%filter_mesh      , &
               u_lon   => dstate%u_lon           , & ! in
               v_lat   => dstate%v_lat           , & ! in
               mfx_lon => block%aux%mfx_lon      , & ! in
               mfy_lat => block%aux%mfy_lat      , & ! in
               mfz_lev => block%aux%mfz_lev      , & ! in
               dmg     => dstate%dmg             , & ! in
               pt      => dstate%pt              , & ! in
               ptfx    => block%adv_batch_pt%qmfx, & ! out
               ptfy    => block%adv_batch_pt%qmfy, & ! out
               ptfz    => block%adv_batch_pt%qmfz, & ! out
               dptdt   => dtend%dptdt            )   ! out
    call block%adv_batch_pt%set_wind(u=u_lon, v=v_lat, mfx=mfx_lon, mfy=mfy_lat, mfz=mfz_lev, m=dmg, dt=dt)
    call swift_prepare(block%adv_batch_pt, dt)
    call adv_calc_tracer_hflx(block%adv_batch_pt, pt, ptfx, ptfy, dt)
    call div_operator(ptfx, ptfy, dptdt)
    call adv_fill_vhalo(pt, no_negvals=.true.)
    call adv_calc_tracer_vflx(block%adv_batch_pt, pt, ptfz, dt)
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds, mesh%full_jde
        do i = mesh%full_ids, mesh%full_ide
          dptdt%d(i,j,k) = -dptdt%d(i,j,k) - (ptfz%d(i,j,k+1) - ptfz%d(i,j,k))
        end do
      end do
    end do
    end associate

    call perf_stop('calc_grad_ptf')

  end subroutine calc_grad_ptf

  subroutine calc_dmgsdt(block, dstate, dtend)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(in) :: dstate
    type(dtend_type), intent(inout) :: dtend

    integer i, j, k

    call perf_start('calc_dmgsdt')

    associate (mesh   => block%mesh   , &
               dmf    => block%aux%dmf, & ! in
               dmgsdt => dtend%dmgsdt )   ! out
    dmgsdt%d = 0
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds, mesh%full_jde
        do i = mesh%full_ids, mesh%full_ide
          dmgsdt%d(i,j) = dmgsdt%d(i,j) - dmf%d(i,j,k)
        end do
      end do
    end do
    end associate

    call perf_stop('calc_dmgsdt')

  end subroutine calc_dmgsdt

  subroutine calc_wedudlev_wedvdlev(block, dstate, dtend, dt)

    type(block_type), intent(inout) :: block
    type(dstate_type), intent(in) :: dstate
    type(dtend_type), intent(inout) :: dtend
    real(r8), intent(in) :: dt

    real(r8) tmp
    integer i, j, k

    ! Follow SB81 vertical advection discretization.

    call wait_halo(block%aux%mfz_lev)

    call perf_start('calc_wedudlev_wedvdlev')

    associate (mesh        => block%mesh           , &
               u           => dstate%u_lon         , & ! in
               v           => dstate%v_lat         , & ! in
               dmg_lon     => block%aux%dmg_lon    , & ! in
               dmg_lat     => block%aux%dmg_lat    , & ! in
               mfz_lev     => block%aux%mfz_lev    , & ! in
               mfz_lev_lon => block%aux%mfz_lev_lon, & ! out
               mfz_lev_lat => block%aux%mfz_lev_lat, & ! out
               dudt        => dtend%dudt           , & ! out
               dvdt        => dtend%dvdt           )   ! out
    call interp_run(mfz_lev, mfz_lev_lon)
    call interp_run(mfz_lev, mfz_lev_lat)
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%full_jds_no_pole, mesh%full_jde_no_pole
        do i = mesh%half_ids, mesh%half_ide
          tmp = -(                                                 &
            mfz_lev_lon%d(i,j,k+1) * (u%d(i,j,k+1) - u%d(i,j,k)) - &
            mfz_lev_lon%d(i,j,k  ) * (u%d(i,j,k-1) - u%d(i,j,k))   &
          ) / dmg_lon%d(i,j,k) / 2.0_r8
          dudt%d(i,j,k) = dudt%d(i,j,k) + tmp
#ifdef OUTPUT_H1_DTEND
          dtend%dudt_wedudeta%d(i,j,k) = tmp
#endif
        end do
      end do
    end do
    do k = mesh%full_kds, mesh%full_kde
      do j = mesh%half_jds, mesh%half_jde
        do i = mesh%full_ids, mesh%full_ide
          tmp = -(                                                 &
            mfz_lev_lat%d(i,j,k+1) * (v%d(i,j,k+1) - v%d(i,j,k)) - &
            mfz_lev_lat%d(i,j,k  ) * (v%d(i,j,k-1) - v%d(i,j,k))   &
          ) / dmg_lat%d(i,j,k) / 2.0_r8
          dvdt%d(i,j,k) = dvdt%d(i,j,k) + tmp
#ifdef OUTPUT_H1_DTEND
          dtend%dvdt_wedvdeta%d(i,j,k) = tmp
#endif
        end do
      end do
    end do
    end associate

    call perf_stop('calc_wedudlev_wedvdlev')

  end subroutine calc_wedudlev_wedvdlev

end module operators_mod
